o 



Methods of robustness analysis for Boolean models 
of gene control networks * 

Madalena Chaves, Eduardo D. Sontag and Rika Albert^ 

<D 
O 
O 

; Abstract 

As a discrete approach to genetic regulatory networks, Boolean models provide an essential quali- 
tative description of the structure of interactions among genes and proteins. Boolean models generally 
assume only two possible states (expressed or not expressed) for each gene or protein in the network as 
well as a high level of synchronization among the various regulatory processes. 

In this paper, we discuss and compare two possible methods of adapting qualitative models to incor- 
porate the continuous-time character of regulatory networks. The first method consists of introducing 
asynchronous updates in the Boolean model. In the second method, we adopt the approach introduced 
by L. Glass to obtain a set of piecewise linear differential equations which continuously describe the 
states of each gene or protein in the network, 
i-^p . We apply both methods to a particular example: a Boolean model of the segment polarity gene 

QT 1 network of Drosophila melanogaster. We analyze the dynamics of the model, and provide a theoretical 

characterization of the model's gene pattern prediction as a function of the timescales of the various 
processes. 

> 

^ ■ 1 Introduction 

in 

Genes and gene products interact on several levels. At the genomic level, transcription factors can activate or 
inhibit the transcription of genes to give mRNAs. Since these transcription factors are themselves products of 
genes, the ultimate effect is that genes regulate each other's expression as part of gene regulatory networks. 
Similarly, proteins can participate in diverse post-translational interactions that lead to modified protein 
functions or to formation of protein complexes that have new roles; the totality of these processes is called 
a protein-protein interaction network. In many cases different levels of interactions are integrated - for 
example, when the presence of an external signal triggers a cascade of interactions that involves biochemical 
reactions, protein-protein interactions and transcriptional regulation. 

During the last decade, genomics, transcriptomics and proteomics have produced an incredible quantity 
of molecular interaction data, contributing to genome-scale maps of protein interaction networks Q |2] [5J 
and transcriptional regulatory networks [4]. Network analysis of these maps revealed intriguing topological 
similarities 0. At the same time, it has been increasingly realized that cellular interaction maps only 
represent a network of possibilities, and not all edges are present and active in vivo in a given condition or in 
a given cellular location j2j[6). Therefore, only an integration of (time-dependent) interaction and activity 
information will be able to give the correct dynamical picture of a cellular network. 
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For many biological networks, and in particular genetic control or regulatory networks, detailed infor- 
mation on the kinetic rates of protein-protein or protein-DNA interactions is rarely available. However, for 
many biological systems, evidence shows that regulatory relationships can be sigmoidal and be well approx- 
imated by step functions. In this case, Boolean models, where every variable has only two states (ON/OFF), 
and the dynamics is given by a set of logical rules, are frequently appropriate descriptions of the network 
of interactions among genes and proteins. Examples include models of genetic networks in the fruit fly 
Drosophila melanogaster 0|U and the flowering plant Arabidopsis thaliana l9l ITOl . 

While Boolean models introduce biologically unrealistic time constraints (typically, such models use 
synchronous updates, which inherently assume that the various biological processes have the same dura- 
tion), they still provide significant qualitative information on the underlying structure of the network. On 
the other hand, continuous models certainly have a more realistic time description of a biological system. 
But, in the absence of information on the kinetic rates, continuous models include many unknown parame- 
ters, and analysis of the system involves exploring the (often large) state space of parameters. An important 
(continuous) model for Drosophila melanogaster segment polarity genes was first developed in [11], where 
a thorough investigation of the parameter space showed that the system is very robust with respect to varia- 
tions in the kinetic constants. Both this continuous model and the discrete model [ 8 1 agree in their overall 
conclusions regarding the robustness of the segment polarity gene network. 

In this paper, we propose two new approaches to the analysis of Boolean models, which combine discrete 
logical rules and structure with more realistic assumptions regarding the relative timescales of the genetic 
processes. The first method introduces asynchronous updates in the Boolean model; since update times 
are randomly chosen, the model is now stochastic (see also [12]). The second method associates to the 
discrete variables a set of continuous variables, whose dynamics is given by a piecewise linear system of 
differential equations, thus introducing a simple "hybrid" model in the manner first suggested by L. Glass 
and collaborators in fl~3l fl4l IT5I . These methods allow us to naturally probe the system with respect to 
perturbations in the time dynamics to analyze its performance. Both methods uncover the robustness of the 
segment polarity gene model [8], and its ability to correctly predict the final gene expression pattern. 

Section |2 summarizes the Boolean model [8], to be analyzed with our methods. The asynchronous 
updates and the Glass-type methods are described and applied in Sections|3]and|4l respectively. In Section|5j 
we discuss the effects of perturbations in initial gene expression, and the development of mutant phenotypes. 
Finally, in Section |6] we show that under a timescale separation between posttranslational and translational 
processes in the model [8], we are able to analytically and exactly compute how frequently the model will 
predict the correct gene pattern. 

2 The segment polarity gene network in Drosophila 

We will apply our analysis to a Boolean model of the interactions among the Drosophila melanogaster seg- 
ment polarity genes. This gene network represents the last step in the hierarchical cascade of gene families 
initiating the segmented body of the fruit fly. While genes in preceding stages of development act transiently, 
the segment polarity genes are expressed throughout the life of the fly. The best characterized segment po- 
larity genes include engrailed {en), wingless (wg), hedgehog (hh), patched (ptc), cubitus interruptus (ci) and 
sloppy paired (sip), coding for the corresponding proteins, which we will represent by capital letters (EN, 
WG, HH, PTC, CI and SLP). Two additional proteins, resulting from transformations of the protein CI, also 
play important roles: CI may be converted into a transcriptional activator, CIA, or may be cleaved to form a 
transcriptional repressor CIR. 

The expression pattern of the Drosophila segment polarity genes (see Table |2]i is maintained almost un- 
modified for three hours, during which time the embryo is divided into 14 parasegments. Each of these 
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Table 1 : Regulatory functions governing the states of segment polarity gene products in the model. Each 
node is labeled by its biochemical symbol and subscripts signify cell number. 



wgi wg i (k+l) = {CIAi (k) and SLPi (k) and not CIRi(k) ) 
or [wgi(k) and (CM<(fc) or SLP^k)) and not CIR^k)] 
Wd WG l {k + l)=wg i {k) 

erii etii(k + 1) = (WGi-i(k) or WG i+1 (k)) and not SLPi(k) 

ENi ENi{k + 1) = em{k) 

hhi hhi{k + 1) = ENi(k) and not CIRi(k) 

HHi HHi(k + 1) = hhi(k) 

ptc { ptc i {k + l) = CIAi (k) and not EN{(k) and not CIRi (k) 
PTd PTd{k + 1) = ptCi(k) or {PTd{k) and not HH^k) 

and not HHi + \(k)) 
ci{ cii(k + 1) = not ENi(k) 
CI, Ch{k + 1) = cii(k) 

CIAi CIAi(k + 1) = CIi{k) and [not PTQ{k) or HH^i(k) 

or HHi + i{k) or hhi_\{k) or hhi + i(k)] 
CIRi CIRi{k + 1) = C7i(fc) and PTCi(k) and not HH^i(k) 

and not HHi + i{k) and not hhi-\(k) and not hhi + \(k) 



parasegments is composed of about 4 cells, delimited by furrows positioned between the wg and en - 
expressing cells [ 16 1. 

The Boolean model that we will study was introduced and developed by one of us in JSJ. (Further 
robustness analysis was also developed in El.) In this model, a parasegment of four cells is considered: 
the variables are the expression levels of the segment polarity genes and proteins (listed above) in each of the 
four cells (the total number of nodes in the network is thus 4x13 = 52.). The expression level of each gene 
or protein is assumed to be either (OFF) or 1 (ON). The model successfully describes the transition from 
the initial expression pattern Q to a final pattern two or three developmental stages later, when the embryo 
has been clearly divided into parasegments (see first entry of Table |2j. As discussed in [8], the evolution of 
these gene expression patterns is well described by a set of logical rules, which are depicted in Tabled 

We adopt the notation "wg\" or "wg 1 (k)" to represent the state of wingless mRNA in the first cell of 
the parasegment at time k. Similar notations apply for other mRNAs and proteins. Periodic boundary 
conditions are assumed, meaning that: node^+i = node\ and node\-\ = node^. The wild type initial 
pattern corresponds to: 



with the remaining nodes zero. 

2.1 Steady states of the Boolean model 

A complete analysis of the steady states is found in [8|. Table |2] summarizes these results, indicating the 
expressed nodes in each of the six steady-states. We note that three of the four main steady states agree 
perfectly with experimentally observed states corresponding to wild type, to ptc knockout mutant (broad 
striped) and to en, wg or hh knockout mutant (non-segmented) embryonic patterns ( 1171 fT8l : see (8| for 
more references). 



Node Boolean updating function (synchronous algorithm) 




wgl = 1, en\ = 1, hh\ 



1) Pt c 2,3,4 ~ 1> c *2,3,4 



1 



(1) 
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Table 2: Complete characterization of the model's steady states. 



Steady state 


Expressed nodes 


wild type 


wg A , WG A , en x , EN X , hh x , HH X , 




ptc 2i , PTC2,3,4, "2,3,4, 




C/2,3,4, CIA 2A , CIR3 


broad stripes 


w §3,4' WGs,i, eni )2 , EN\ >2 , hh 1>2 , HH 1>2 , 




ptc 3i , PTC 3A , c/ 3j4 , C/ 3) 4, CIA 3A 


no segmentation 


"1,2,3,4? C/1^,3,4, -fTCl,2,3,4> CZRi 2, 3 ,4 


wild type variant 


wg 4 , WG 4 , eni, ENi, hh x , HH X , 




ptC 2 4 , PTCi t 2,3,<i> "2,3,4, 




C/2,3,4, C/A 2 , 4 , C//? 3 


ectopic 


wg 3 , WG3, en 2 , EN 2 , hh 2 , HH 2 , 




ptc 13 , PTC lj3)4 , cii j3) 4, 




C/1,3,4, C/A13, CIR4 


ectopic variant 


wg 3 , WG 3 , en 2 , EN 2 , hh 2 , HH 2 , 




ptc l 3 , PrCi )2)3 ,4, "1,3,4, 




C/1,3,4, C/A13, CIR4 



2.2 The regulatory function of the sloppy paired gene 

The rule for SLP protein in Table [J summarizes in a simple way the experimental observations on the ex- 
pression and regulatory activity of the sloppy paired gene in the segment polarity network [19]. A more 
detailed rule for the sloppy paired expression pattern can be created to incorporate recent evidence of en- 
grailed protein inhibiting sip transcription [20|. However, inhibition by engrailed accounts only partially 
for the experimentally observed restriction of sip to the posterior half of the parasegment. Thus we need to 
invoke an additional regulatory effect, which we denote by RX. RX probably represents a combination of 
regulation by the pair-rules responsible for the establishment of sip, namely runt, opa and Factor X [21 1 and 
of sip autoregulation. 

Therefore, SLP expression in Tablefjcan be replaced by the following set of equations: 



RXi(k + l) 



JO, if i € {1,2} 

\ 1, if ie{3,4} 

(2) 

slpi (k + 1) = RXi (fc) and not ENi (k) 
SLPi(k + l) = slpi(k) . 

The sloppy paired initial conditions would then be: 

slp° 3A = 1, SLP° 3A = 1. (3) 

This generalization of the segment polarity network model introduces additional steady states, such 
as a two-stripe en pattern characterized by slp 4 = SLP 4 = 1, wg 4 = WG4 = 1, en\ j3 = ENi t3 = 1, 
hh\ :3 = HH\ 3 = 1, ptc 2 4 = PTC 2 4 = 1, and c/2,4 = C/2,4 = CIA 2 4 = 1 (this pattern was also found 
in [ 8 1 as a result of sip mutation). This expression pattern is non- viable since it has two en — wg borders and 
would lead to an ectopic parasegment structure. On the other hand, it is not difficult to see that, starting from 
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conditions Q and (J3}, none of the "new" steady states are reachable, since these initial conditions imply 
that neither sip nor SLP can change their expression at any time. Indeed we have the following result: 

Lemma 2.1 Consider the extended model of Table ^together with ©. Assume that initial conditions are 
given by © and ©. Then slp^t) = sip® and SLPi(t) = SLP° for all times . 

A sketch of proof is as follows. Note first that slp l 2 (t) = SLP\^{t) = follows from RX\^{t) = for 
all t. Next, observe that: 

slp 4 (tx) = => ENi(t 2 ) = 1=> en 4 (t 3 ) = 1 =>- SLP 4 (t 4 ) = =>- slp 4 (t 5 ) = 0, 

where t\ > t% > t% > t 4 > £5 > 0. Thus, in order for slp 4 to become zero at any time, it had to be so 
at some previous instant (t§ < t±). If slp 4 (0) = 1, and if T > is the first instant such that slp 4 (T) = 
then we have a contradiction. Hence slp 4 (t) = 1 for all times. (Note that this argument is independent 
of the order in which the nodes are updated.) Similar arguments show that slp 3 (t) = slp 3 (0) = 1 and 
SLP 3A (t) = 1, for all t. 

Thus the extended model leads to the same results as assuming a constant SLP pattern. For this reason, 
and lacking more specific biological evidence on the regulation of sloppy paired, our present analysis is 
focused on the simpler biologically relevant model of Tabled 

3 Asynchronous algorithms 

In general, for a network of N gene products (denoted x\, . . ., xn), the dynamics of a Boolean model is 
typically studied by simultaneously updating the state of all the nodes in the network, according to 

x! +1 = Fi{Xtxl...,Xk), i = l,...,N 

where Fi is the regulating function for mRNA or protein X, L . An underlying hypothesis is the existence 
of perfect synchronization among the various regulatory processes. However, it is well known that the 
timescales of transcription, translation, and degradation processes can vary widely from gene to gene and 
can be anywhere from minutes to hours. 

In analogy with task coordination and data communication procedures in the context of parallel compu- 
tation systems \22\, we have previously developed several methods that introduce different timescales for 
the different regulatory processes within the network [ \2\. These include algorithms that randomly choose 
the order in which the nodes are updated (this random order algorithm is summarized below in 13. II and in 
Section|6ll, and a totally asynchronous algorithm, where the next updating times for each node are randomly 
chosen at each instant. 

We now introduce a more intuitive asynchronous algorithm, where each node is updated according to its 
own specific time unit. The time units for the nodes are randomly chosen from a uniform distribution in an 
interval [1 — e, 1 + e], where e G (0, 1). The updating times of i-th node are then pre-specified as: 

Tl,T?,...,T?,... fcGN, 

with 

T k+l = T k + 1 . = klh k£N _ (4) 

For instance, / y WG4 < j wgi means that wingless protein in the fourth cell is translated at a faster rate (shorter 
time intervals) then wingless mRNA is produced. At any given time t, the next node(s) to be updated is(are) 
i such that T* = minj/{Tj > t}, for some k. The variables X{ are updated according to: 

X t (T l k ) = F l (X 1 (T? i ),...,X N (T* i )), (5) 
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where rj| defines the most recent instant when node j was updated, that is 

r^ = max{T/: T*<I*}. (6) 

By ordering all the time sequences {T k : i = 1 . . . N, k = 1, 2, . . .}, into a single nondecreasing sequence, 
say {t\,t2, ■ ■ •}, the asynchronous model can also be written in the form 



Xi(tk+l) 



Fi(Xi(t k ),X 2 (tk), ■ ■ ■ ,X N (t k )), ififc = hi for some i 
Xi(tk), otherwise. 



It is clear that the steady states of this model must satisfy Xi(tk+i) = -Xj(ifc) = Fi(X(tk)), and therefore 
are the same as those of the synchronous boolean model in Table |2] 

Note that the case e = reduces to the synchronous model, where every node is updated simultaneously 
(7i = 1), at the same time instants: T k = k, for all i = 1, . . . , N. This algorithm allows great variability in 
each process' duration, exploring the gene expression patterns due to all possible combinations of individual 
timescales. Implementation of this asynchronous algorithm shows that, if started from the initial wild type 
state 0, any of the steady states of the model (Tabled may occur with a certain probability. The probability 
of occurrence of each pattern depends on the range over which the individual time units 7j are allowed to 
vary (see Fig. |2j. For e = 0, the wild type steady state is attained with probability 100% (corresponding 
to the synchronous Boolean model). As e increases to 0.01 (resp. 0.1) this value decreases to 60% (resp. 
44%). However, further increase in e (hence larger time intervals) unexpectedly leads to an increase in 
the occurrence of the wild type state, up to 51% for s = 0.9. Other final states observed are the broad- 
striped pattern (25% — 38%) observed in heat-shock experiments and ptc mutants ifTTl and the pattern with 
no segmentation (12% — 15%) observed in en, hh or wg mutants lfl~8l . the latter two corresponding to 
embryonic lethal phenotypes 1171 . Each of the other three steady states occurs with frequencies less than 
5%. (These values were obtained from 10000 numerical experiments.) 

A possible extension of this algorithm would be to consider a discrete model with a finite number of 
logical levels describing ON, OFF as well as other intermediate steps of the system [7 1. This would involve 
decisions about the number of intermediate steps, their values, and development of new transition rules. 
Instead, in this paper we will focus on a "hybrid" model, that takes into account the continuous nature of the 
biological processes, while still using Boolean rules to describe ON/OFF transitions (Section©. 



3.1 Random order algorithm 

For comparison purposes, we briefly summarize an alternative asynchronous algorithm, which guarantees 
that every node is updated exactly once during each unit time interval [ 12 1. A random order of updates for 
the N nodes is generated as a permutation (j) k of {1, . . . N}. This permutation is randomly chosen out of a 
uniform distribution over the set of all ./V! possible permutations, at the beginning of the time unit k. The 
updating times for each node are now written as 

= N{k-i) + (j) k {i), ken, 

so that 4> k {j) < 4> k (i) implies T k < T k , and node j is updated before node % at the k-th iteration. The 
results of this algorithm are qualitatively similar to those of the asynchronous algorithm (Tabled- 



4 Glass-type networks 

The asynchronous algorithm defined by ©, (|5} and © allows the introduction of distinct timescales for 
each regulation process in a Boolean model. We next propose an alternative method which provides a 
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Table 3: The frequencies of the six steady states observed with the three different methods when starting 
from the wild type initial condition. 



Steady state 


Asynchronous 


Random order 


Glass-type 


pattern 


algorithm 


algorithm 


model 


wild type 


44-51% 


56% 


89-100% 


broad stripes 


25-38% 


24% 


0-6% 


no segmentation 


12-15% 


15% 


0-3% 


wild type variant 


4-5.6% 


4.2% 


0-1% 


ectopic 


0.4-1% 


0.98% 


0% 


ectopic variant 


0.1-0.5% 


0.68% 


0% 



bridge between discrete and continuous approaches, resulting in a more realistic model, but without the 
necessity of specifying any kinetic or binding parameters (which are typically unknown). In this method, 
the gene and protein levels are represented as continuous variables, and their time evolution is described 
by differential equations, but the interactions among nodes are still modeled by Boolean functions H3|[T4l 
El IIU- Glass [ 13 1 introduced a class of piecewise linear differential equations that combine logical rules 
for the synthesis of gene products with linear (free) decay by describing each node with two variables, one 
discrete and one continuous. For simplicity of notation, in what follows we will let X,; denote the continuous 
variable associated with node i, its discrete variable Xj, and the discrete variable's Boolean rule by Fj. The 
Glass-type model is then 

dX- — 

—^ = -X i + F i (X 1 ,X 2 ,...,X N ), i = l,...,N. (7) 
at 

At each instant t, the discrete variable X{ is defined as a function of the continuous variable according to a 
threshold value: 

where 9 G (0, 1). The discrete variables X{ represent the ON and OFF levels of the nodes in the Boolean 
model. The underlying assumption in this Glass model is that the decay rate and activation threshold of 
each gene product is identical. Since the initial condition for the piecewise linear system Q is also Q 
(i.e., X(0) = X(0)) and Fj G {0, 1}, it is easy to see that solutions of Q evolve in the hypercube [0, 1]^. 
Under these conditions, the limiting values "0" and "1" of the continuous variable Xi represent, respectively, 
"absence of species i" and "maximal concentration of species i" - thus we can view the X as dimensionless 
variables, scaled to attain their maximal values at 1. The continuous dynamics is translated into a Boolean 
ON/OFF response, according to 9: as soon as X increases above 9, species i is considered to be in the ON 
state; otherwise it remains in the OFF state (see also [23]). Thus the parameter 9 defines the fraction of 
"maximal concentration" necessary for a protein or mRNA to effectively perform its biological function. 
This method allows us to study the continuous evolution of the genetic network simply by specifying 9, 
the fraction of maximal concentration that is effective as ON level, avoiding the need to specify any kinetic 
parameters. Below and in Section 16.41 we will see that system Q exhibits distinct dynamics in the two 
regions 9 < 1/2 and 9 > 1/2. It is easy to see that the steady states of the piecewise linear equations Q are 
still those of the Boolean model, since: 

^1 = 0^ X=Fi(X u X 2 ,...,X N ), i = l,...,N, 
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independently of 9. Applying this method to the Drosophila segment polarity gene network, we find an 
exact convergence to the wild type steady state when started from the wild type initial condition (Fig. 0, 
independently of the ON/OFF threshold value 9. This result supports the Boolean model as a suitable 
description of the underlying network of gene interactions. 

4.1 Introducing distinct timescales 

The assumption of equal decay rates and activation thresholds for all nodes is an oversimplification similar 
to that made in synchronous Boolean models. However, for further robustness analysis, one may intro- 
duce different timescales for the different processes, by scaling the time units in each differential equation 
according to: 

dX- — 

-^ = a i (-X i + F i (X 1 ,X 2 ,...,X N )), (9) 

with 04 > e for some fixed e > 0, i = 1, . . . , N. Each Xi is a discrete variable defined as before (note that 
steady states of this new system are still those of the Boolean model). 1 

This method represents a continuous equivalent to the asynchronous algorithm described in Section |3] 
Here, the (inverse) scaling factors a^ 1 may be viewed as half-lives of mRNA or proteins. These may be 
directly compared to the individual time units 7, as follows. Using Euler's method to discretize system Q 
obtains: 

%{t + At) = %{t) + aiAt(-%(t) + Fi(X{t))). 

Now notice that choosing the integrating time interval to be such that ctjAi = 1 recovers the discrete 
asynchronous algorithm with specific time units 

7i = At = a' 1 . 

For comparison to the discrete algorithm, we choose both the scale factors a" 1 and the time units 7$ ran- 
domly from a uniform distribution in intervals of the form [1 — e, 1 + e], e € (0, 1). The numerical experi- 
ments are reported in Fig. 13 where the threshold value 9 ® was set to 1/2. We again observed that, starting 
from wild type initial conditions Q, all steady states may occur with a certain frequency (see TableO. But, 
in contrast to the asynchronous Boolean model, the wild type pattern occurs with frequencies that decrease 
monotonically with e, down to 89% for e = 0.9 (Fig. |2]>. The next more frequently achieved patterns are the 
broad stripes (Fig. @), with probability 6% for e = 0.9, the no segmentation (Fig.15}, with probability 3%, 
and the wild type variant, with probability 1%. 

The three methods we have described (Table |3j produce qualitatively compatible results, in the sense 
that the wild type pattern is always the most frequently occurring steady state, followed by the broad stripes, 
no segmentation, and wild type variant patterns. 

4.2 Fraction of maximal concentration that defines an ON state 

The piecewise linear system © follows the threshold (JSJl to decide whether a given node is ON or OFF. 
While this value 9 did not affect the dynamics of the system in the case a\ = a 2 = • • • = ocjy, it plays 

'We will analyze the behaviors of trajectories of systems of the form 0, assuming that trajectories are well-defined. Since the 
right-hand sides of equations of these type are discontinuous, it is very difficult to give general existence and uniqueness theorems 
for solutions of inital-value problems. One must impose additional assumptions, insuring that only a finite number of switches can 
take place on any finite time interval, and often tools from the theory of differential inclusions must be applied, see for instance 1 24 
for more discussion. See also 1 23 1. 
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a significant role in the general case. In Figure |3j it is immediate to see that the effect of 9 depends on 
the length of the interval allowed for the timescales. Indeed, there is a marked difference between narrow 
(e < 0.5) and wide (e > 0.6) intervals. For narrow intervals, numerical experiments indicate that, starting 
from initial condition (Q, system © converges to wild type steady state with probability around 90% or 
more; whereas, for wider intervals, this probability may decrease down to 68% (e = 0.9, 9 = 0.9). 

On the other hand, the threshold value also divides the dynamics into two regions. For 9 < 0.5, the 
probability that initial condition O leads to the wild type steady state is above 90%, independently of e. For 
9 < 0.5 it is less probable to reach the non segmented (around 1%) than the broad striped pattern (around 
10%). For higher 9 > 0.6, the probability of reaching the wild type steady state (from initial condition Q) 
decreases very significantly, and in addition, it becomes more problable to reach the non segmented (around 
20%) than the broad striped pattern (around 7%). 

In our Glass-type model, 9 represents the fraction of maximal concentration above which an mRNA 
or protein is considered ON, or biologically effective. Our results indicate that quite small fractions of the 
maximal concentration can (and should) be interpreted as sufficient amounts for an mRNA or protein to be 
in the ON state. Namely, when the concentration of mRNAs or proteins has increased to a fraction up to half 
its maximum possible value, it is already present in a sufficient amount to perform its function. This follows 
from observation of Figure|3j < 9 < 0.5 leads to a good (realistic) performance of the model, with 90% 
convergence to wild type pattern. Our results also indicate that a higher threshold is, on the contrary, not 
very realistic. For instance, setting 9 > 0.6 leads to a fairly high incidence on the mutant patterns. But by 
letting 9 > 0.6 one is assuming that an mRNA or protein is not present in a sufficient amount to perform its 
function until it reaches at least sixty per cent of its maximal value. Typically such high thresholds are not 
observed: indeed, in ifTTTl (supplementary material) where continuous dynamics are also transformed into an 
ON/OFF response, a threshold of 10% is used, and considered very reasonable. In our experiments, unless 
otherwise indicated, we have used 9 = 0.5. 

To further generalize system (|9j, it is natural to allow each variable to respond to its own threshold and 
consider distinct 9i values. We will address this problem in Section l6~3l and see that the system's dynamics 
is preserved in each 9 region. In fact, our simulations with distinct 9i recover many of the theoretical results 
we obtained for a species-independent 9 (Section l6~4t . 

5 Pre-patterning errors and knockout mutant situations 

In |8l El we identified some sufficient or necessary initial conditions for obtaining the wild type steady 
state (minimal pre-patterns). We now analyze how mutant patterns arise from gene "knockout" experiments 
or delay in establishing the initial pre-pattern. 

Gene knockout experiments consist of completely supressing the expression of a given gene in all cells. 
In our models this is equivalent to setting the corresponding mRNA permanently zero in all equations. Thus 
a wingless knockout can be analyzed by setting wg i = for all i = 1,2, 3, 4, and in every equation of the 
Boolean rules in Tabled The steady states of the resulting system can now be computed from: 

X* = F*(X*), (10) 

where the vector X* and functions F* include all but the knockout variables. For example, for the wingless 
knockout, we drop the functions F wg . and for the other mRNAs and proteins we have 

WGi = F; c . =0, i = l,2,3,4 

node = F n * odc = F node , for nodes en, hh, ci,ptc, EN, HH, PTC, CI, CIA, CIR, 
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Table 4: The steady states corresponding to gene knockouts in the segment polarity network model, calcu- 
lated according to (11 Ob , Here, for the ci knockout, the "wild type" state is interpreted as the wild type pattern 
in all but the ci, CI, CIA,CIR mRNA/proteins. 



Knockouts 


Mutant steady states 


wg, en, hh 


no segmentation 


ptc 


broad stripes 




wild type, 


ci 


broad stripes , 




ectopic 



because wg appears explicitly only in the rules F WGi . So, it follows immediately that WGi = 0, and the 
remaining nodes are then also easy to compute. It is easy to check (see also [8]) that knockouts of wg, en, 
hh or ptc exhibit only one steady state, while knockouts of ci exhibit three steady states (summarized in 
Table© in both asynchronous and Glass-type models. 

From another point of view, one may consider a delay in the establishment of the pre-pattern (that is, 
the full initial condition Q). If expression of a given gene is delayed, does the system recover, and how 
soon? To answer this question, we simulated a delay in expression of gene X, by setting the corresponding 
discrete variable X(t) = in all cells, for all t < T May . We then varied T dclay between and 7 time units, 
and measured the frequency of occurrence of each steady state, both for the asynchronous and Glass-type 
models. (In the latter we set the concentration threshold equal to 1/2 for all nodes.) 

The results are shown in Figures |6] and We can see that, for short T dcky , both models recover their 
original frequencies of occurrence of each steady state, while for long T delay both models converge to the 
corresponding mutant steady state. A curious exception is the case of ci, where long delays in its initial 
expression do not significantly change the probability that the wild type steady state is achieved (and even 
slightly increase it in the asynchronous model). This agrees with the conclusion of [ 8 ] that cubitus inter- 
rupts knockout coupled with an otherwise wild-type initial condition converges to a state close to the wild 
type steady state. Remarkably, delays in both cubitus interruptus transcription factors (CIA, CIR) have a 
lesser effect than an imbalance in their expression (see fT2l ). This leads us to predict that, during the pre- 
segmentation stage of embryo development, the cubitus interruptus proteins' expression is the last to be 
established. 

Another noteworthy observation is the fact that small delays in wg expression have much more drastic 
effects on the system than the same delay in en or hh. This phenomenon reflects the one-way signaling 
cascade starting with expression of wg, which induces en, which in turn induces hh. We see that a total 
disruption in the system is caused by a delay of only 3 time units in wg or en expression, which cause the 
system to fail to reproduce the wild type pattern, and settle into a non-segmented pattern. However, if only 
hh is delayed, the system is disrupted only after a delay of 5 time units. In other words, recovery of the 
system back to the "good" developmental process is more probable in the event of a hh expression delay, 
than a wg or en expression delay. 

6 Robustness of the model under timescale separation 

In the previous algorithms, the space of all possible timescales for protein/mRNA regulatory processes was 
explored, with no assumptions on the characteristic duration of translational or post-translational processes. 
As a consequence, the robustness analysis shows that the model diverges from the wild type pattern very 
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often, with the biologically inviable states occurring with a noticeable frequency. However, it is also well 
known that post-translational processes such as protein conformational changes or complex formation, usu- 
ally have shorter durations than uanscription, translation or mRNA decay. This fact justifies the introduction 
of a distinct timescale separation among processes, by choosing to update proteins first and mRNAs later. 

6.1 Timescale separation in the random order algorithm 

Timescale separation is straightforwardly implemented in the random order algorithm presented in Sec- 
tion |3j at the k-th updating step we generate two random permutations, <p Pml and 4>^ lRNA , within the set of 
proteins and mRNAs, respectively. Then the N nodes are updated in the order given by 

4> k = (4L, 4> k mRNA )- 

This method again shows that the Boolean model is very robust, in the sense that when started from the wild 
type initial condition, the wild type pattern occurs with a frequency of 87.5% and only one other steady 
state is observed, the broad striped pattern, with a frequency of 12.5%. Furthermore, these frequencies are 
exact as we show in lfT2l . where we also completely characterize the model resulting from incorporation of 
a protein/mRNA timescale separation into the random order algorithm. We show that the wild type state is 
in fact an attractor for the system, while the pathway to the broad stripes state may exhibit oscillatory cycles. 
We summarize this and other results in the next theorem, stated without proof, and refer to [ 12 1 for more 
details. 

Theorem 1 In the random order algorithm with timescale separation, let wg 3 = 0, ptc 3 = 1, hh® 4 = and 
c/g = 1 (as satisfied by the initial pattern Q). Then system diverges from the wild type pattern if and only 
if the permutation (f) 1 satisfies the following sequence among the proteins CI, CIA, CIR and PTC: 

CIR 3 Ch CIA 3 PTC 3 , 

Ch CIR 3 C1A 3 PTC 3 , (11) 

CI 3 CIA 3 CIR 3 PTC 3 , 

while the other proteins may appear in any of the remaining slots. I 

Thus we can compute the exact probability with which the random order algorithm (with timescale 
separation) leads to either the wild type or broad stripes pattern: the latter is simply the fraction of sequences 
of the form {n}, and equals 12.5% El. 

6.2 Timescale separation in the Glass-type and asynchronous algorithms 

As shown in Section |4j the (discrete) asynchronous algorithm and the (piecewiese continuous) Glass-type 
system provide equivalent representations of a gene expression network. Indeed, the "specific time units" 
7i and the inverse "scaling factors" a" 1 , both represent the rate of dynamical evolution of each individual 
node. For these two models, we implement time separation among processes by using two non-overlapping 
intervals for the scaling factors: 

7~\aj e A mRNA , if Xi E {wg,en,hh,ptc,ci} 

7 t ~\ai £ A Pmn if Xt e {WG, EN, HH, PTC, CI, CIA, CIR}, 

with, for instance, A mRNA = [0.2, 0.6] and A Pm = [1.4, 1.8]. Under these conditions, choosing the factors 
from a uniform distribution in these intervals, numerical experiments indicate that the two methods respond 
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Table 5: Probabilities of convergence to a given steady state, under the separation of timescales assumption. 
These values are theoretically exact for both the random order algorithm and Glass-type model. 



Steady state 


Random order 


Asynchronous 


Glass-type 


pattern 


algorithm 


algorithm 


model 


wild type 


87.5% 


93.7-100% 


100% 


broad stripes 


12.5% 


0-6.3% 


0% 


other 


0% 


0% 


0% 



in mostly similar ways, with only two patterns occurring at steady state when the systems start from (wild 
type) initial condition ©. The two possible steady states are the wild type and broad stripes patterns. 

For the asynchronous algorithm, the probabilities of convergence to each of the steady states clearly 
depend on the distance between the two intervals: convergence to wild type is 93% in the case where 
intervals A mRNA , A Prot are consecutive, and up to 100% for the case 2a < b, a G A mRNA , b G A Pml . 

For the Glass-type model, two cases can be distinguished: 9 < 0.5 and 9 > 0.5. For 9 < 0.5, numerical 
simulations show that the model reaches wild type pattern with probability near 100%, even when there 
is some overlap between A mRNA and A Pmt (see Fig. |8j. In fact, we next theoretically prove that the wild 
type pattern is indeed the unique possible steady state of the hybrid system © and initial condition ©, as 
indicated by the simulations, when there is a suitable distance between the intervals, and a lower bound on 9 
(TheoremEJ Section WM . For 9 > 0.5, we have found no condition that guarantees convergence to the wild 
type steady state, and indeed numerical simulations show that, even for large interval separation, the system 
may converge to one of the mutant patterns. 

A comparison of Theorems^and|2]emphasizes differences and similarities between discrete and contin- 
uous models: intuitively, the single discrete event described by Theorem^cannot take place in a continuous 
model. Therefore wg 3 remains "0" (OFF) for all times, ruling out the possibility that the broad stripes pat- 
tern is reached. Indeed, Theorem ^establishes that (in the discrete case, with the random order algorithm) 
divergence from wild type pattern occurs if and only if wg\ = 1. This fact involves a jump in wg 3 from "0" 
to "1" at precisely the first iteration. On the other hand, in the Glass-type model, the continuous variable wg 3 
cannot instantaneously jump from "0" to "1". Since the discrete ON/OFF levels are defined by a threshold 
on wg 3 , there will necessarily be a smoothing effect on any transition between "0" and "1". This is what 
happens in case (a) of Theorem 13 

The second (sufficient) condition of Theorem |2] guarantees convergence to the wild type steady state 
for all < 6 < 0.5 (while condition (a) is only for 0.382 < 9 < 0.5), but assumes that a PTC3 > a c , A . 
This is an analog to Theorem [l] if a PTCi > a ci3 , then (starting from PTCs(0) = Cls(0) = and assuming 
Fptc 3 = Fa s = 1) PTC% increases faster than C/3, implying that PTC3 becomes ON faster than C/3. Such 
response prevents the events listed in Theorem [0 which would lead to a mutant state. Thus, both discrete 
and piecewise linear model predict that the sequence of PTC, CI expression in the third cell is one of the 
fundamental pieces in establishing the correct development of embryo segmentation. 

It is also worth pointing out that the three methods provide qualitatively similar results under timescales 
separation, all predicting that only wild type and broad stripes patterns to occur, the latter with considerably 
smaller frequency. 

6.3 Distinct concentration thresholds for ON state 

A natural question arising in the analysis of equation © concerns the dynamics of the system under more 
general concentration thresholds. We have seen that 9 (even when equal for all species) plays an important 
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role in establishing basins of attraction to each of the steady states of the model. 
We have, in particular, identified three distinct regions of behavior: 

Region 1 : < 9 < 
Region 2 : ^fi < 9 < \ 
Region 3 : \ < 9 < 1 

When 9 is equal for all nodes, the probability that the system evolves into the wild type pattern is above 90% 
in regions 1 and 2, but may be as low as 68% in region 3. Furthermore, for region 2, we theoretically prove 
that that probability is exactly 100% under the timescale separation assumption. 
We now associate to each node a specific 9i, so that (J8jl is modified to: 

* w \ i, ut) > e, . 

To test the performance of the system and compare it to previous results, we considered two timescale 
situations: G [0.5,1.5] for all i, or the timescale separation A mRNA = [0.2,0.6], A Pm , = [1.4,1.8]. In 
each case, we randomly assigned values to 9i from uniform distributions in the intervals (0, 1), (0, 0.5) and 
(0.4, 0.5). (Note that 0.4 is close to (3 — \/5)/2.) The numerical results with varying 9j extend and confirm 
our previous observations for 9i = 9, i = 1, . . . , N. 

Table[6]summarizes the results for all combinations of 9i and a» regions. The most general case, allowing 
a large degree of freedom in both timescales and concentration thresholds, indicates the vulnerability of the 
network, with a very high incidence on mutant patterns (#j G (0, 1), ai G [0.5, 1.5]). Again we see that 
there is a marked difference in the 9{ regions below or above 0.5. Comparing the reasonable results obtained 
for 9i < 0.5 with the bad performance for 0.5 < 9{ < 1, we conclude that the optimal ON concentratin for 
proteins or mRNA is below 50% of maximal concentration. 

Restricting 9i even further to one of the conditions given in Theorem l2l (Section 16.41 conditions devel- 
oped for the case 0j = 9, i = 1, . . . , N) dramatically increases the probability that the system develops in 
the correct way. 

These simulations suggest, moreover, that some of the theoretical results obtained for each of the three 
regions may extend to the case of distinct 9{. Indeed, note that, under timescale separation, when 9{ are 
chosen from region 2, convergence to wild type steady state is 100% - compare to part (a) of Theorem|2] 

These simulations further confirm the role of a PTC3 and a c , A in the segmentation network. Requiring 
a PTC3 > a> ci3 decreases the probability of formation of the broad stripes pattern (any timescales), but doesn't 
influence the probability of a non segmented embryo. The latter mutant is prevented only by a complete 
separation of timescales of the regulatory processes. 

6.4 Glass-type model provides exact convergence to wild type pattern 

In this section we will require that the intervals A mRNA and A Pm , do not overlap, by satisfying the following 
assumption: 

For all a G A mRNA and b G A Pm , : < 2a < b. (13) 

A second assumption is that the effective maximal concentration is equal for all nodes and satisfies 9 < 1/2, 
which is equivalent to: 

1 1 

In <ln- (14) 

l-9~9 y 
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Table 6: Probabilities of convergence to a given steady state, with distinct concentration thresholds 0, and 
distinct timescales on, in the Glass-type model. (Probabilities computed out of 1000 simulations for each 
case.) 



Steady state 


(0,1) 


(0.5,1) 


(0,0.5] 


(0,0.5] 


[0.4,0.5] 


0.5 


0i 


pattern 
















wild type 


45.6% 


57.1% 


84.1% 


90% 


92.6% 


94.2% 


A — 

■"■mRNA — 


broad stripes 


27.8% 


15.1% 


12% 


6.2% 


7.3% 


4.5% 


Ap m , = 


no segmentation 


24.4% 


25.8% 


0.9% 


0.9% 


0.05% 


1.3% 


[0.5,1.5] 


wild type, variant 


2.1% 


1.9% 


2.9% 


2.9% 


0% 


0% 




wild type 


74.1% 


52.7% 


96.6% 


97.1% 


100% 


100% 




broad stripes 


10.8% 


3.3% 


3.3% 


2.8% 


0% 


0% 


[0.2,0.6] 


no segmentation 


14.1% 


43.9% 


0% 


0% 


0% 


0% 


Ap rBI = 


wild type, variant 


1.0% 


0% 


0.1% 


0.1% 


0% 


0% 


[1.4,1.8] 



Theorem 2 Consider system © with initial condition Q. Assume that the scaling factors on satisfy ill 3b . 
Assume also that one of the following conditions holds: 

(a) 9 < 1/2 and (1 - 9) 2 < 9 or equivalently 0.382 « (3 - y/S)/2 < 9 < 1/2; 

(b) 9 < 1/2 and a PTC3 > a ci3 ; 
then wg 3 (t) = for all t. 

This shows that the steady state representing the broad stripes pattern cannot ever be reached in system © 
from the initial condition Q, when 6 < 1/2 and either of the extra conditions holds. 

Theorem 3 Consider system © with initial condition Q. Assume that the scaling factors aii satisfy (fl3t . 
and that 9 < 1/2. Then wg A (t) = 1 and PTC x {t) = for all t. 

This shows that the steady states represented by the no segmentation, wild type variant or the two ectopic 
patterns also cannot ever be reached in system © from the initial condition Q. From Theorems |2]and[3]we 
conclude that, under the timescale separation assumption, the Glass-type model © can only converge to the 
wild type pattern, when starting from the initial condition Q, and for appropriate 6 values (Tabled. 

We first summarize some useful observations. Let X denote any of the nodes in the network, and a its 
time rate. Since equations © are either of the form dX/dt = a(— X + 1) or dX/dt = — aX, their solutions 
are continous functions, piecewise combinations of: 



x\t) = 1 - (1 -X 1 (t ))e- a( ^* o) (15) 
X°(t) = X°(t )e- a{t - to) (16) 



(resp. X (i)) is monotonically increasing (resp. decreasing). In addition, note that discrete variables 
X can only switch between and 1 at those instants when X(i swiKh ) = 9, that is: 

tLi = k+ ^i^M> ( ,7> 

a 1 — 

= *o + -m^ (18) 
a 9 
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From the initial conditions, together with the constant values of SLPi = 1,2,3, 4), we can immediately 
conclude: 



wg lj2 (t) = WG 1>2 (t) = 
eh 3A (t) =EN 3A (t) = 0, 
hh 3A (t) =HH 3A (t) = 0, 



(19) 



(20) 



for all £ > 0. Then, because ci 3A (0) = 1 and F d3 4 = not EN 3 a, 



ci 3A {t) = 1 and C7 3j4 (£) = 1 



(21) 



Lemma 6.1 Let < £o < t 3 < £i and < £ 2 < t 3 . Define 5 = In / max^...^ «i- Assume CIA 3 (t) = 
for £ G (t 2 ,t 3 ), and wg 3 (t) = for t £ [0, t 3 ). Then 

(a) wg 3 (t) = for t G [0,t 3 + <ty 

(b) WG 3 (i) = for t e [0, t 3 + 5); 

(c) en 2 (i) = ^2(t) = for t £ [0, t 3 + 5); 

(d) M 2 (t) = HH 2 (t) = for t G [0, t 3 + 5). 
Assume further that PTC 3 (t) = 1 for i e (t , Then 

(e) PTC 3 (t) = 1 for all i G (t , t 3 + <5). 

(f) C/A 3 (t) = for all t G (i 2 , t 3 + 

Proof: Part (a) follows directly from the fact that F wg3 (t) = on [0, t 3 ), and from (fT7l . 

To prove parts (b), (c), and (d), first note that initial conditions together with wg 3 (t) = for t G [0, t 3 ) 



for t G [0, £3]. Then, from equations dT5l to (fT8t we conclude that the corresponding discrete variables 
cannot switch from to 1 during an interval of the form [0, £ 3 + ^- In j^q)- Taking the largest common 
interval yields the desired results. 

To prove parts (e) and (f), assume also that PTC 3 (t) = 1 for £ G (£o,£i)- From d20t and part (d), it 
follows that function F PTC . A does not switch in the interval (£0, £ 3 + 5) and in fact PTC 3 {t) = 1 for all £ in 
this interval. This, together with d20l) and part (d) yield F CM3 (£) = for (£o,£ 3 + 5), so that CIA 3 cannot 
increase in this interval and the discrete level satisfies CIA 3 (t) = for all £ G (£2, £ 3 + 5), as we wanted to 
show. I 

Corollary 6.2 Let < £ < £ 3 < £1 and < £ 2 < £ 3 . If PTC 3 (t) = 1 for £ G (£ ,£i), CIA 3 (t) = for 
t G (£ 2 , £3), and wg 3 (t) = for £ G [0, £ 3 ), then wg 3 (t) = for all £. 

Proof: Applying Lemma l6~Tl we conclude that, given any k > 0: 



imply 



WG 3 (t) = 0, a 2 (£) = HV 2 (£) = 0, 
hh 2 (t) = HH 2 (t) = 0, 



CIA 3 (t) = 0, for £ G (£ 2 ,£ 3 + kS) 
wg 3 (t) = 0, for £ G [0, £ 3 + kS) 
PTC 3 {t) = 1 for £ G (£ ,£ 3 + k8) 
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imply 



CIA 3 (t) = 0, for te(t 2 ,t 3 + (k + 1)5) 
wg 3 {t)=0, for te [0,t 3 + {k + l)5) 
PTC 3 (t) = 1 for t E (t ,t 3 + (k + 1)5). 

Since 5 is finite, we conclude by induction on k that wg 3 (t) = for all t. 
Proof of Theorem^ The rule for CM 3 may be simplified to (by d20t ) 

F CIAa = CI 3 and [notPrC 3 or hh 2 or HH 2 \. 

From equation dm . we have that 



CI 3 (t) = 1, for all t > — In -. (22) 

a c/3 1-0 



On the other hand, since ptc s (0) = 1, by continuity of solutions ptc 3 (t) = 1 for all t < ^— In i. This 
implies that the Patched protein satisfies 

PTC 3 (t) = 1 - e ~ aprc 3 *, < t < — In - 



ptC3 



and therefore 



f 0, < t < ln-r^j 

™ 3(t) = ^-r» <( ^mj. (23 > 

By assumption, Q/» rC3 > a pK3 and also In < In a, defining a nonempty interval where PTC3 is expressed. 
Now let t c = In and t p = — — In jKt. CIA 3 {t) starts at zero and must remain so while CI 3 = 0, 
so that 

CZA 3 (t) =0 for < t < t c . 

In the case t c > t p , letting to = t p , t± = In ^, t 2 = 0, and £3 = t c in Corollary 16.21 obtains wg 3 (t) = 
for all t. This proves item (b) of the theorem, and part of (a). 

To finish the proof of item (a), we assume that (1 — 9) 2 < 9 and must now consider the case t c < t p . 
Then 



CIA 3 (t) 



0, < t < t c 

tp 

_1 

Following equation dT7l with to = tc an d CIA 3 (to) = 0, CIA3 might become expressed at time t c <t a < t„: 



l-e _acM 3 (*-*<=), t c <t<t p 
clA 3 (t p )e-—A^ P ), tp <t<^-\nl, 



1 , 1 

t a = t c -\ In ■ 



«C/A 3 1 — 9 

but it would then become zero again at (equation dT8b with to = t p ) 

1 CIA 3 (t p ) 
tb = t p -\ In 



a 



cias 
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Finally, we show that, even if CIA^(t) = 1 for t G (t a , %), wg 3 cannot become expressed in this interval. In 
this interval, wg 3 evolves according to wg 3 (t) = 1 — e~ Q, ' s 3 (* - *«) ) and wg 3 can switch to 1 at time 

1 , 1 

t w = t a -\ In • 



a wg3 1-9 

We will show that t w > tf,, so wg 3 (t) = in the interval [0, Writing 

CM 3 (t p ) _ CIA 3 (t p )l-6 



In ^-^ = In 



1 

CM 3 (t p ) , 1 



= In + In 

1-8 

In — — + In 



where we have used CIA 3 (t p ) < 1 and the assumption on 9: < jKj. Therefore 

2 , 1 

h < t v H m ■ 



up -r i a 

&CIA 3 1 — 
1,1 1,1 

< In r + In < U 

a„„„ 1 — a 



CIA 3 17 

where we have used the timescale separation assumption ( fT3t . Letting to = t p , t<i = 0, and ti = t 3 = 
min{tfe, aQ* In |} in the Corollary, obtains wg 3 (t) = for all t. I 
We will next show that if wg 4 (t) = 1 in a given interval [0, T), then in fact wg 4 (t) remains expressed for 
a longer time, up to T + 5, with 5 > 0. This is mainly due to assumption (fT3l . which says that mRNAs take 
longer than proteins to update their discrete values, because they have longer half-lives: a^ NA > a~} r This 
allows the initial signal "wg 4 = 1" to travel down the network, sequencially affecting the wingless protein, 
engrailed, hedgehog and CIA, and feed back into wingless allowing wg 4 to remain expressed for a further 
time interval. 

Lemma 6.3 Let T > In I and define 

— a„ g4 V 

a wc A i 

In ^ — '-. (24) 



a WG4 u 

If wg A {t) = 1 for < t < T, then 

(a) WG 4 (t) = 1 for t G In -L,T + 5); 

(b) eni{t) = 1 for t G [0,T + <5); 

(c) SVi(t) = 1 - e - Q£ V for i G [0, T + <5), and HVi(t) = 1 for i n y + $); 

(d) cii(t) = 0, Ch{t) = 0, CMi(t) = 0, and CJRxif) = for f G [0, T + 5); 

(e) = 1, fort G [0,T + <5); 

(f) C/A 4 (t) = 1, for t G In ^ + ^- In T + 5), and C//? 4 (t) = 0, for t G [0, T + 5); 

(g) wg 4 (i) = lfortG [0,T + $). 
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Proof: Let T > In ^, and assume that wg 4 (i) = 1 for < t < T. To prove part (a), note that 

WG±{t) is of the form <TT3T> (with to = 0, and M7 4 (0) = 0) and the corresponding discrete variable is 
WGAt) = 1, for t € In jK, T). Moreover, suppose that wgAt) = for t > T, then 

WG 4 (i) = (1 - e - QM/G 4 T )e" a,VG 4( t - T ), t > T. 

But WG 4 remains 1 until the switching threshold is attained, that is up to time 

1 (1 _ e - a wa 4 T\ 
T+ ln^ 1 

Ot W G 4 V 

> T+ In i '- 



a 



H'6'4 



= T + 5. 

Thus we conclude that WG&(t) = 1 in the desired interval. 

To prove part (b), observe that F mi (t) = WG 4 (t) for all t, from ill 9b . and recall that eni(0) = 1. From 
part (a), F enl (t) = 1 for t E (^~~ ^ n i~h?> + ^ n tne °th er hand, en\ can only switch from 1 to at 
t = In p which is larger than a~g 4 In So, in fact, en\(t) = 1 for all < t < T + 5. 

Part (c) follows immediately by integration of the EN± equation. 

To prove part (d), first recall F dl = not ENi and the initial conditions ci\ (0) = = CI\ (0) = C1A\ (0) = 
Cfl?i(0). Therefore ci\(t) increases up to t = In and then decreases in a^ N \ In < t < T + 5. 

Now note that the discrete variable ci±(t) remains in the whole interval [0, T + 5). This is because ci\ never 
reaches the threshold: this would be attained at some t > In but, since oQ 1 In > In jK,, 
the function ci\ starts decreasing before it could reach the value 0. Finally, from the rules of the Cubitus 
proteins it is immediate to see that Ch{t) = CIA\{t) = CIRi(t) = for t £ [0, T + 5). 

To prove part (e), recall that F m = EN\ and not CIR\. From part (a), it follows that F hhl (t) = 
in the interval [0, acj^ In j^) and F m (t) = 1 in the interval (a~^ In jz^,T + 5). Since Mi(0) = 1, 
hh\(t) decreases in the interval [0, In jK*) but increases in (a^J In , T + 5). The discrete value is 
hh\{t) = 1 in the whole interval, since hhi(t) remains above the threshold. (The justification is similar to 
the case of ci\ (t) in part (d).) 

To prove part (f), note that part (e) and then the use of d2TT >. allows us to simplify F CIAi : 



F CIA4 (t) = Ch{t) and hh x (t) = 1, t G ( In- ^T + 5). 

a c , A 1-9 



Thus 



C/A 4 (t) 



0, < t < -±- In 



—ariA. t In -r— s , . 

1-e 4 V aa * **J, -^ln^ <t<T + 5, 



and CYA 4 (t) = 1 for t E [^- In + a In ^— ^, T + 5). Observe that this interval is indeed nonempty, 
by assumption O- Finally, F clRi (t) = CI 4 (t) and notMi(i) = 0, and hence CIR 4 (t) = for t G [0, T + 5). 
To prove part (g), we note that (from part (f)) 

F wg4 (t) = l, te(—ln-L- + J_ln-l-,r + «5), 
a c/4 1-0 a CM4 1-0 
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implying that vvg 4 (t) increases in this interval. On the other hand, we know that vvg 4 (i) > 9 and wg A (t) = 1 
up to at least t = In I > In H — In j\ . This shows that in fact we At) = 1 for all 

r a„ (4 6 aa 4 1-9 a C iA 4 1-8 * 4V ' 

t€[0,T + S). I 

Proof of Theorem^ Since wg 4 (0) = 1, from equations (117b . ill 8b . we know that the earliest possible 
switching time from 1 to is In ^. Applying Lemma lo31 with T = a~* In ^ establishes that wg 4 (i) = 1 
for t E [0, T + <5), with (5 given by <24l . Next, applying LemmalOwith T = a~ x In ^ + k5, k <E N, shows 
that wg A (t) = 1 for i G [0, T + (fc + 1)<5). Since 5 is finite, we can conclude by induction that wg A (t) = 1 
for all t > 0. 

To prove that PTCi(t) = 0, note that CIA\{t) = (Lemma IO with T = +oo) implies ptc x (t) = 0. 
Since PTC\(0) = and PTC\ cannot become expressed unless ptc 1 is first expressed, the desired result 
follows. I 



7 Conclusions 

We discussed two alternative methods for modeling gene expression networks: purely discrete Boolean 
methods and piecewise linear differential systems, which combine continuous degradation with discrete 
synthesis. For both methods we introduced new techniques for a deeper analysis of the networks with 
respect to perturbations in the timescales of the system. For the piecewise linear system we also studied the 
effect of the ON concentration thresholds. 

We find that unrestricted variability in the duration of the diverse processes present in the network 
may lead to significant deviations from experimentally observed results, thus suggesting the fragility of the 
developmental process under severe perturbations (the asynchronous algorithm fails to predict the correct 
pattern with probability 50%, and the Glass-type system fails with probability at least 10%). 

Another set of numerical experiments introduces a separation between timescales of post-translational 
and transcription/translation processes, and in practice "updates mRNAs later than proteins". In this context, 
the piecewise linear Glass-type system indicates a remarkable robustness of the Boolean model in predicting 
the final gene expression pattern. Indeed, we provide a theoretical proof that the Glass-type system always 
correctly generates the wild type development (i.e. the convergence to the wild type steady state when started 
from the wild type initial state), under the separation of timescales assumption, and appropriate OFF/ON 
thresholds. The asynchronous model's predictions depend very much on the degree of separation between 
timescales. As the intervals A Pmt and A mRNA become closer, the asynchronous algorithm increasingly fails to 
generate the wild type pattern. This leads us to conclude that a strong separation between timescales of post- 
translational and transcription/translation processes is necessary for establishing the regular gene expression 
pattern in the segment polarity network. From analysis of the piecewise linear model we conclude that the 
fraction of maximal concentration above which a protein or mRNA is effectively ON needs to be quite small, 
below 50%. Higher concentration thresholds may disrupt the development process. 

The comparison between discrete and continuous models shows clearly that sudden transitions may 
happen in discrete systems and lead to a false result; such sudden transitions are smoothed out in continuous 
models which prevent generation of false results (see also [25 1). Nevertheless, we conclude that both models 
agree in predicting the fundamental sequences of gene expression that irreversibly lead to a deviation in the 
development towards a mutant state. 

By combining continuous-time techniques with discrete events, we can with great generality explore 
and sample the space of all possible timescales as well as of effective ON levels. Moreover, as information 
about the mRNA/protein lifetimes, decay rates or activation thresholds becomes available, it can be straight- 
forwardly incorporated by fixing the corresponding inverse scaling factor ol 1 . The hybrid model retains 
the ease of Boolean models in determining the steady states corresponding to gene knockouts and perturbed 
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initial conditions. It is straightforward to calculate which mutant patterns result from each gene knockout. 
Here we also studied the the effect of perturbations on the prepattern, by simulating delay in initial expres- 
sion of each gene. We find that the system is vulnerable to large delays (larger than two time units) in 
expression of any gene - except for ci - and, in such delayed conditions, the mutant state characteristic to 
that gene knockout is generated. For low order delays, the system typically recovers and proceeds through 
the correct wild type development. 

The Glass-type system with time separation, as a model of the segment polarity gene network, reflects 
the conclusion of von Dassow et al. that the topology of the network is more important than the fine-tuning 
of the kinetic parameters ifTTl . since its results are robust for a large region of parameter (scaling factor, ac- 
tivation threshold) space. Due to its underlying Boolean structure, the model also intrinsically incorporates 
the recent finding of Ingolia that parameter sets need to satisfy certain constraints - that ensure the bista- 
bility of certain genes- to lead to correct solutions [ 26 ]. Taken together, the results of the synchronous [8|, 
asynchronous |12| Boolean and hybrid models convincingly demonstrate the Boolean models' capability 
for effectively describing the basic structure and functioning of gene control networks when detailed kinetic 
information is unavailable. 
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Figure 1 : The solution to the system of piecewise linear equations (Q (dashed lines). Each column represents 
one cell, and each of the 13 rectangles represents the continuous-time variable for proteins or mRNAs (as 
labeled at left). In each rectangle, the y-axis ranges from to 2 units. The nodes for which the trajectories 
converge to 1 (middle of the rectangle) are exactly those expressed in the wild type steady state. (The time 
units are arbitrary.) 
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Figure 2: Probability of occurrence of the three most frequent patterns: wild type (WT), broad stripes (BS), 
and no segmentation (NS), under variable range of timescales. Dashed lines/squares represent asynchronous 
algorithm results, while solid lines/circles represent Glass-type model results (out of 1000 runs). Results 
were obtained with 9 = 0.5. 
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Figure 3: Effect of effective ON concentration, 6, on the probability of occurrence of the three most frequent 
patterns: wild type (WT), broad stripes (BS), and no segmentation (NS). The results are for the Glass-type 
model with aj 1 randomly chosen in an interval [1 — e, 1 + e]. Solid lines represent the case e = 0.5, dotted 
lines represent the case e = 0.1 and dash-dotted lines represent the case e = 0.9. (results out of 1000 runs). 
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Figure 4: A solution to the system of piecewise linear equations © (dashed lines), in an example where the 
steady state corresponds to the broad stripes pattern. In each rectangle, the y-axis ranges from to 2 units. 
Notice that wingless is expressed in two adjacent cells, as opposed to the wild type pattern, where wingless 
is expressed in only one cell (similarly for engrailed and hedgehog). 
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Figure 5: A solution to the system of piecewise linear equations © (dashed lines), in an example where 
the steady state corresponds to the no segmentation pattern. In each rectangle, the y-axis ranges from to 
2 units. Notice that wingless, engrailed and hedgehog are not expressed in any cell, thus no segments are 
visibly detected in the embryo. 
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Figure 6: Effect of initially delayed expression in occurrence of steady state patterns. Left: delay in wingless 
expression. Middle: delay in engrailed expression. Right: delay in hedgehog expression, (results out of 
1000 runs). 



27 




Figure 7: Effect of initially delayed expression in occurrence of steady state patterns. Left: delay in patched 
expression. Right: delay in cubitus interruptus expression, (results out of 1000 runs). 
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Figure 8: Probability of occurrence of the three most frequent patterns: wild type (WT), broad stripes (BS), 
and no segmentation (NS), with separation of timescales. Dashed lines/squares represent asynchronous 
algorithm results, while solid lines/circles represent Glass-type model results. The x-axis represents the 
level of separation, computed by min{6 G A^,}/ max{a 6 A mRNA } (out of 1000 runs). Results were 
obtained with 6 = 0.5. 
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